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Abstract 

An accurate computational method is presented to determine the mass distribution in a rotating thin- 
disk galaxy from given rotation curve by applying Newtonian dynamics for an axisymmetrically rotating 
thin disk of finite size with or without a central spherical bulge. The governing integral equation for mass 
distribution, resulting from the balance between the Newtonian gravitational force and centrifugal force 
due to rotation at every point on the disk, is transformed via a boundary-element method into a linear 
algebra matrix equation that can be solved numerically for any given rotation curve. The mathematical 
formulation of the thin-disk model can easily be extended to including a central spherical bulge. To 
illustrate the effectiveness of this computational method, mass distributions in several mature spiral 
galaxies consistent with various types of measured rotation curves are determined without the need of 
fictitious rotation velocity outside the "cut-off" radius. When a central spherical bulge is present, the 
total galactic mass increases only slightly but the mass distribution in the galaxy is altered in such a 
way that the periphery mass density is reduced while more mass appears toward the galactic center. 
By extending the computational domain beyond the galactic edge, we can determine rotation velocity 
outside the cut-off radius which appears to continuously decrease and gradually approach the Keplerian 
rotation velocity out over twice the cut-off radius. In examining the circular orbit stability, the galaxies 
with flat or increasing rotation velocities with radius seem to be more stable than those with decreasing 
rotation velocities especially in the region near the galactic edge. 

keywords: galaxy: disk — galaxies: general — galaxies: kinematics and dynamics — galaxies: struc- 
ture — methods: numerical and analytical 

1 Introduction 

Without direct means for accurate measurement, mass distribution in galaxies — gravitationally bound 
assemblies of (10 5 -10 12 ) stars — can only be inferred from the observable information according to the 
known physical laws. In astronomy, the observable information is usually carried by electromagnetic 
radiation — the light — emitted from the visible objects. The light can be analyzed to provide informa- 
tion about the emitting objects such as their material constituents, surface temperature, distance, moving 
velocity, etc. Observations have shown that many (mature spiral) galaxies share a common structure 
with the vis ible matter distributed in a flat thin disk, rotating about their center of mass in nearly circular 
orbits (e.g. jBinnev & Tremainelll987l) . The speed of circular motion of objects in galaxies can be deter- 
mined from the measured Doppler shift of light and its plot against the radial distance from the galactic 
center is called the rotation curve or circular speed curve. Using the measured rotation curve has been 
considered as the mo st reliable means for deriving mass distribution in thin-disk galaxies dToomrel 19631 
Sofue & R ubin 200 j). 
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Another independent means for estimating mass distribution is based on luminosity measurements 
of galactic surface-brightness profile by assuming a given (usually constant) mass-to-light ratio, the va- 
lidity of which seems to be rather debatable especially when it comes to quantitative estimation of mass 
distribution. In fact, descrepancies often arise between the observed rotation curves of galaxies and that 
predicted from mass distributions following the surface-brightness profile based on Newtonian dynam- 
ics, leading to the so -called "galaxy rotation probl em" as still haunting the astrophysical community to 
the present day (e.g jFreeman & McNamaral2006h . If not being relied on for quantitatively determining 
the mass distributi on, the typically exponential deca ying profile of observed surface brightness in many 
galaxies (Freeman 119701 Binnev & Tremaine 1 19871) suggests the general structure of decreasing (sur- 
face) mass density with the radial distance from the galactic center, which appears quite consistent with 
that derived f rom the measured ro tation curves according to Newtonian dynamics for rotating thin-disk 
galaxies (e.g jFeng & Galloll201 at least in a qualitative sense. 

Given a measured rotation curve, to derived mass distribution in a thin-disk galaxy also requires the 
physical laws that can make the connection between the kinematic behavior and locations of matter. 
For galactic dynamics, th e best known, well-establis hed physical laws are Newton's laws of motion and 
Newton's law of gravity (Bin nev & Tremainelll987l) . Thus, we focus herewith on the mass distribution 
in rotating thin-disk galaxies derived from measured rotation curves according to Newtonian dynamics. 

Although theoretically well-established, the actual computations of Newtonian dynamics when ap- 
plied to thin-disk galaxies appeared to be much more difficult than that for a gravitational system with 
spherical symmetry such as the solar system. Serious efforts were ma de for integrating the Pois - 
son equation with mass sources di s tributed on a disk, as summarized by iBinnev & Tremaind (119871) . 
Bratek . Jalocha & Kutscheral(l2008h . lFeng & Ga5old201 ll) . among others. The solution directly obtained 
from such efforts is usually the (Newtonian) gravitational potential which can yield the gravitational 
force by taking its gradient. In an axisymmetric disk rotating at steady state, the gravitational force (the 
radial gradient of gravitational potential) is expected to equate to the centrifugal force due to rotation at 
every point. 

Unlike the spherically symmetric mass distribution that generates the gravitational force at a given 
radial position only depending upon the amount of mass within that radius, the gravitational force due to 
a thin-disk mass distribution can be influenced by matters both inside and outside that radius. Thus, the 
mass distribution in a thin-disk galaxy cannot be determined simply by applying Keplerian dynamics that 
relates the mass within a radial position to the rotation speed at that radius. In principle, the rotation speed 
at a radial position is mathematically related to mass distribution in the entire disk of the galaxy. The 
fact that the brightness in disk galaxies typically decreases exponentially with radial distance indicates 
a practical limitation of the rotation curve measurements: the detectable signal must terminate at a 
finite radial position — the so-called "cut-off radius". All measured rotation curves terminate at their 
cut-off radii, although sometimes the cut-off radii may move further out with new signal detection and 
processing technology development. 

Among several possible approaches of solution, using Bessel functions has been the method of 
choice for many authors (e.g., iToomrell 1 9631 iFreemanl 1 970t iNordsieckll 1973b ICuddefordl 1993 ; Conwav 



20001: 1 Jalocha. Brateck & Kutscherall2008l:lBratek, Jalocha & Kutscherall2008l) probably due to the con 
venience in theoretical derivations. The mathematical formulations with Bessel functions typically con- 
tain integrals extending to infinity, which has become the major practical difficulty when working with 
the rotation curves that always terminate at finite cut-off radii. The part of rotation velocity outside 
the cut-off radius must be constructed based on various assumptions wi thout much measurable infor- 
mation, to complete the mathemat i cal formulation (as discussed by 



Nordsieck 1973; Bosma 1978 



Jalocha, Brateck & Kutscheral 2008: Bratek, Jalocha & Kutschera 2008). 

To avoid the need of the part of fictitious rotation curve outside the cut-off radious, an integral equa- 
tion for a rotating thin-disk galaxy with its edge coinciding with the cut-off radius of rotation curve can 
be formulated according to Newtonian dynamics, c onsisting of Green's function in terms of the complete 
elliptic integrals of the first kind and second kind (IFeng & GalloluOl lb . With appropriate mathematical 
treatments, the apparent numerical diffic ulties associated wit h singularities in elliptic integrals can be 
completely removed (as demonstrated by Feng & Gallo 201 1 , when carefully evaluating the mathemat- 
ical limit). To enable dealing with arbitrary forms of rotation curves and mass density distributions, the 
boundary element method for solving integral equations is adopted here using compactly supported basis 
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functions instead of that extending to infinity like Bessel functions. Hence the finite physical problem 
domain for a disk with edge ending at finite radius can be conveniently considered by solving a linear 

algebra matri x problem. 

Following lFeng & Gallol d201 ll) . in the present work we nondimensionalize the governing equations 
such that a dimensionless parameter, which we call the "galactic rotation number", appears in the force 
balance (or centrifugal-equilibrium) equation, representing the ratio of centrifugal force and gravitational 
force. Together with a constraint equation for mass conservation, the value of this galactic rotation num- 
ber can be determined as part of the numerical solution. The value of the galactic rotation number can be 
used for determining the total gala ctic mass in the disk from measured galactic (cut-off) radius and char- 
acteristic rotation velocity. While lFeng & Gallol (1201 lb focused mainly on illustrating the computational 
method with a few somewhat idealized rotation curves, here we apply this method to in-depth analysis of 
the realistic rotation curves available in the open literature (e.g., http://www.ioa.s.u-tokyo. ac.jp/~sofue, 



de Blok et alj|2008l , etc.) We also extend our method to including spherical central core and bulge, to 



further applications such as for determining rotation velocity beyond the cut-off radius, and so on so 
forth. 



2 Mathematical Problem Description 

For convenience of mathematical treatment, a rotating galaxy is represented by a self-gravitating con- 
tinuum of axisymme trie ally distributed mass in a circular disk with an edge at finite radius R g (beyond 
which we expect mass density to diminish precipitously to the inter-galactic level having inconsequential 
gravitational effect on the galactic disk dynamics). Without loss of generality, we consider the thin disk 
having a uniform thickness (h) with a variable mass density (p) as a function of radial coordinate (r). In 
the situation of thin disk, the vertical distribution of mass (in the z-direction) is expected to contribute 
inconsequential dynamical effect especially as the disk thickness becomes infinitesmal. In mathematical 
terms, the meaningful variable here is actually the surface mass density er(r) = p(r) h. Whether to 
consider the surface mass density a{r) or the bulk mass density p(r) in the mathematical equations is 
just a matter of taste, since they can easily be converted to each other using a constant factor h by our 
definition without substantial difference. Here, we choose to use the bulk density p(r) for its consistency 
with the common physical perception of a thin disk with a nonzero thickness h. 

For steady rotation, there must be a balance between the gravitational force and centrifugal force 
at every point in the galactic disk. If the force density on a test mass at (r, 9 — 0) generated by 
the gravitational attraction due to the summation (or integration) of a distributed mass density p(f) at 
position described by the variables of integration (f, 9) is expressed as an integral over the entire disk, 
with the distance between (r, 9 = 0) and (?, 9) given by (f 2 + r 2 — 2rr cosf)) 1 ' 2 and the vector 
projection given by (f cos 9 — r), the equation of force balance in a rotating thin disk can be written as 
(according to Newton's laws) 



(f cos 9 — r)d9 
o (f 2 + r 2 — 2fr cos 9) 3 / 2 



p(f)hfdf + A 



V(r) 



0, 



(1) 



where all the variables are made dimensionless by measuring lengths (e.g., r, f , h) in units of the out- 
ermost galactic radius R g , disk mass density (p) in units of Mj/Rg with denoting the total mass 
in galactic disk, and rotation velocities [V(r)] in units of the a characteristic galactic rotational velocity 
Vq (usually defined according to the rotation curve of interest). The disk thickness h is assumed to be 
constant and small in comparison with the galactic radius R g . The numerical results for surface mass 
density p(r) h are expected to be insensitive to the exact value of the ratio of h/R g as long as it remains 
small. There is no difference in terms of physical meaning between the notations (r, 9) and (f, 9); but 
mathematically the former denotes the independent variables in the integral equation (fTJ whereas the 
latter the variables of integration. The gravitational force represented as the summation of a series of 
concentric rings is described by the first (double integral) term while the centrifugal force by the second 
term in ((T). 

Nondimensionalizating the force-balance equation yields a dimensionless parameter, which we call 
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the "galactic rotation number" A, as given by 



^ Vp R g 
M d G 



(2) 



where G (= 6.67 x 10~ n [m 3 /(kg s 2 )]) denotes the gravitational constant, R g is the outermost galactic 
radius which may be taken as the cut-off radius of measured rotation curve, and Vb is the characteristic 
velocity. This galactic rotation number A simply indicates the relative importance of centrifugal force 
versus gravitational force. 

Equation (Q]) can either be used to determine the surface mass density p(r) h from a given rotation 
curve V(r) or vice versa. But when both p(r) and A are unknown, another independent equation is 
needed to have a well-posed mathematical problem. In view of the conservation of mass, the total mass 
of the galaxy disk should stay as a constant satisfying the constraint 



2ir I p{f)hfdr = 1. 
lo 



(3) 



This constraint offers an addition equation for determining the value of galactic rotation number A. With 
equations (HJ-OJ, the mass density distribution p(r) in the disk, the galactic rotation number A, and the 
disk galactic mass Md can all be determined from the measured values of V(r), R g , Vb, and h. On the 
other hand, if p{r) and h (or p(r) h) as well as A are given, V(r) can of course be determined from ([T}. 
The integral with respect to 8 in (Q} is known to be equivalent to 



(r cos 8 — r)dO 



E(m) K(m) 



"(f — r) r(r + r) 



(f 2 + r 2 — 2fr cos #) 3 / 2 

where K(m) and E(m) denote the complete elliptic integrals of the first kind and second kind, with 

4rr 



(4) 



Thus, (Q]i can be written in a single-integral form 

~E{m) K(m) 



I 

Jo 



(f + r) 2 



p(f)hfdf 



-AV{rf 



0. 



(5) 



(6) 



which is more suitable for the boundary element type of n umerical implementatio n 



Following a standard boundary element approach (e.g., Sladek & Sladek 1998; Sutradhar. Paulino & Gray 



2008), the governing equations © and (O can be discretized by dividing the one-dimensional problem 
domain [0, 1] into a finite number of line segments called (linear) elements. Each element covers a sub- 
domain confined by two end nodes, e.g., element i corresponds to the subdomain [r^, fj+i], where and 
ry+i are nodal values of r at nodes i and i + respectively. On each element, which is mapped onto a 
unit line segment [0, 1] in the ^-domain (i.e., the computational domain), p is expressed in terms of the 
linear basis functions as 

p(0 = /»i(i-o+Pi+i£, o<e<i, (7) 

where pi and pi + i are nodal values of p at nodes i and i + respectively. Similarly, the radial coordinate 
f on each element is also expressed in terms of the linear basis functions by so-called isoparameteric 
mapping: 

f(0=h(i-0 + n+^, o < e < i. (8) 

If the rotation curve V(r) is given (as from measurements), the N nodal values of pi = p(r^) are deter- 
mined by solving N independent residual equations over N — 1 element obtained from the collocation 
procedure, i.e., 



N-l 

E 

71=1 



E(mi 



(IV 



i 



-AV(r 



= 0,i=l,2,...,JV, 



(9) 
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with 

m«(fl= fr, ( l )ri (10) 

where = /3 n (l — £) + Pn+iC- The value of A can be solved by the addition of the constraint equation 

N-i r i 

CiT 

P^MO^df- 1 = 0. (11) 



2-E 



Thus, we have A + l independent equations for determining A+l unknowns; the mathematical problem 
is well-posed. With appropriate mathematical treatments of the singularities arising from the elliptic 
integrals and boundary conditions at r = and r = 1 (as described in Appendix A), the set of linear 
equations (0 and ( fTTT i for A + 1 unknowns (i.e., N nodal values of pi and A) c an be put in a m atrix 
form and then solved with a standard matrix solver such as by Gauss elimination ( Press et al1ll988l) . 



3 Results 



To obtain numerical solutions, the value of (constant) disk thickness h must be provided; we assume 
h = 0.01 out of many possible choices. For computational efficiency, we distribute more nodes in the 
regions (e.g., near the galactic center and disk edge) where p varies more rapidly. Unless the rotation 
curve has very steep velocity changes that need finer discretization with more elements, the typical 
number of nonuniformly distributed nodes N used in computating most cases is 1001 (corresponding to 
1000 linear elements) with which we found for most cases to be sufficient for obtaining a smooth curve 
of p versus r and discretization-insensitive values of galactic rotation number A. 

The rotation curves available in the open literature (e.g., http://www.ioa.s.u-tokyo. ac.jp/~sofue) are 
typically provided in a tabular form with data points at radial positions often not coinciding with our 
nodal positions. We use the cubic spline interplation method (IPress et al.Ml988D to evaluate our nodal 
values of V(r) from the rotation curve data such that the rotation curve used in our computations is 
guaranteed to smoothly pass through all the measured data points. 



3.1 NGC 4736 



The galaxy NGC 4736 has recently been studied bv lJalocha. Brateck & Kutsch era (2008). for illustrating 
that the baryonic matter distribution can account for the observed rotation curve. Thus, we believe it 
deserves our attention of analysis using our computational method. 

There are several different versions of rotation curve data for NGC 4736 in the literature. Here we 
consider two of them, one is from the websit e of Sofue (e.g., h ttp://www.ioa.s.u-tokyo. ac.jp/~sofue) and 
the other from THINGS measurements (Ide Blok et al.ll2.Q08b . Figure 1 shows the two versions of the 



rotation curves with r measured in units of R g = 10.35 (kpc) (where 1 kpc = 3.086 x 10 19 m), and 
rotation velocity V(r) in units of Vn = 150 (km/s). 

As shown in iFeng & Gallol d201 lb . the value of total galactic mass in the disk can be determined 
according to (O with computed value of A as 



Ma 



vjR 9 

AG 



(12) 



Because the computed value of the galactic rotation number A is 1.9656 for the THINGS rotation 
curve and 1.5908 for that of Sofue, we obtain Ma = 2.756 x 10 10 (solar-mass) (where 1 solar-mass 
= 1.98 8 9 2 x 10 30 kg) when the THINGS rotation curve is used and M d = 3.405 x 10 10 (solar-mass) 
when the Sofue rotation curve is used, for the NGC 4736 galaxy. The value of Md = 3.405 x 10 10 
(solar-mass) agrees well with that computed by IJalocha. Brateck & Kutschera! ( 12008b (i.e., 3.43 x 10 10 
solar-mass) u sing the same rotation curve of Sofue . 

However. IJalocha. Brateck & Kutschera! ( 2008 ) used an iterative spectral method with Bessel func- 
tions which requires the inclusion of rotation curve beyond the cut-off radius extending to infinity. They 
also considered the mass density due to hydrogen HI outside the cut-off radius. With our method, only 
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Figure 1: Profiles of NGC 4736 rotation velocity V(r) and mass density p{r), with the thick line for that of 
THINGS and the thin line for that from Sofue. The computed values of the galactic rotation number A are 
1.9656 and 1.5908 for the case of THINGS and Sofue, respectively. 
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Figure 2: The distributions of mass density p(r) with (thin line) or without (thick line) hydrogen HI extend- 
ing beyond the cut-off radius, and the corresponding rotation curves. Only the case based on the THINGS 
rotation curve is shown here. 



7 



the available data for rotation curve within the cut-off radius is needed. Meanwhile, we assume the mass 
density in the galactic disk diminish at the same cut-off radius to enable a self-consistent consideration 
of the mathematical problem on a finite disk domain. The solution of the axisymmetric mass distribution 
in the galactic disk for a given rotation curve by one-step G auss elimination of the linear algebra matrix 



equation without further iterations (cf . iFeng & 



If desired, the effect of hydrogen HI and the unmeasurable rotation curve outside the cut-off radius 
can be conveniently exami ned in an a posteriori manner. For example, the surface mass density of 



hydrogen HI considered by Jaloch a. Brateck & K utschera (2008) at the cut-off radius is about one solar- 



mass per square pc, which translates to our nondimensional p = R^/(Md h) — 0.3887 (or 0.3146, with 
R g measured in units of pc and Md in units of solar-mass) for the THINGS (or Sofue) rotation curve, 
decreasing one order of magnitude in about 3 (kpc) beyond the cut-off radius at 10.35 (kpc). 

To examine the effect of hydrogen HI outside the cut-off radius, we modify the mass distribution 
starting from a radius n < 1 such that the mass density for r > r\ is described by 

p(r) = p ie -[0-2+(— 1)/0-3] 2 , ri < r < oo , (13) 

where r\ — 0.965 and p\ = p{r{) — 0.388. The profile of mass density distributi on extending outside 



the cu t-off radius (r = 1) described by ( fT3] > approximates well to that considered bv Lfalocha. Brateck & Kutscheral 
(2008) while simplifying the analysis. With the given mass distribution extending beyond r = 1, we can 
corrrespondingly extend the integration beyond r = 1 in © and (0, to calculate rotation velocity be- 
yond the cut-off radius. Figure 2 shows the mass density distribution with (thin line) and without (thick 
line) the hydrogen HI outside the (nondimensional) cut-off radius r = 1, and the corresponding rotation 
curves. The integration result of (0 shows that including hydrogen HI beyond the (nondimensional) 
cut-off radius r = 1 increases the total galactic mass only by ~ 0.5%. Therefore, it is not surprising 
to notice that the original rotation curve in figure 2 is altered insignificantly by the mass density modi- 
fication according to Sl3[ . In fact, the rotation curve beyond r = 1 calculated without the mass density 
modification Sl3[ differs so little from that in figure 2 that it is visually almost indistinguishable when 
plotted together. 

In view of typical uncertainties in rotation curve measurements (as illustrated in figure 1 for two 
different versions of the same galaxy), we expect that any matters (such as hydrogen HI) outside the 
cut-off radius of a galaxy cannot have substantial influence on the gross galactic rotation characteristics, 
because the amount of mass in comparison to that of Md is ususally so insignificant. Thus, whether 
including the hydrogen HI mass outside the cut-off radius of NGC 4736 should have inconsequential 
effect on Newtonian dynamics relating the measured rotation curve to mass distribution in the galactic 
disk. As illustrated here, however, consideration of HI mass beyond r = 1 can be conveniently im- 
plemented as an a posteriori process (without iteratively computing solutions) to evaluate (or, in other 
word, predict) the rotation velocity beyond r = 1 which could not be obained from measurements. Then, 
the needed part of rotation curve beyond the cut-off radius, for using a formulation that requires it (e.g. , 
Nordsieck|[l97l iBosmall 19781; IJalocha. Brateck & Kutscherall2008tlBr"atek. Jalocha & Kutscherdl2008l) 
to determine mass distribution from measured rotation curve, can be provided using our method with 
concrete certainty (namely, without fictitious assumptions). 



3.2 Milky Way, NGC 4945 

The Milky Way, also called the Galaxy, is the galaxy that contains the Sun and the Earth, which is why 
it is of particular interest to astronomy and astrophysics. NGC 4945 is a spiral galaxy that appears quite 
similar to the Milky Way. So, we present results for both of them together here. The rotation curve 
data provided by Sofue (e.g., http://www.ioa.s. u-tokyo.ac.jp/~sofue) suggest nonzero rotation velocity 
at r = 0. According to our continuum treatment of a rotating disk galaxy with Newtonian dynamics 
description of force balance ©, nonzero rotation velocity at r = requires a strongly singular mass 
density to ensure that 

l 

p{f)df — s- oo . (14) 
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Figure 3: Profiles of the Milky Way and NGC 4945 rotation velocity V(r) and mass density p(r), with the 
thick line for that of the Milky Way and the thin line for NGC 4945. The computed values of the galactic 
rotation number A are 1.6365 and 1.6873 for cases of the Milky Way and NGC 4945, respectively. 
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This is because that the kernel of integral in © for any nonzero f has a limit value of zero at r = 0, i.e.. 

' E[m) ir(m)1 - 1 [^(0)-iv(0)] , (15) 



lim 

r->0 



and E(0) = K(0) = tt/2 dAbramowitz & Stegunlll972h . Thus, V(Q) must be zero according to © 



unless (TPiT l is true ( as in the disk oflMesteJ 19631 where p — > 1/r as r — > 0). 



As discussed in Fens & Gallo (1201 II) . the computational method used here can reproduce the result 
of Mestel's disk dMeste]||l963l) for the entire problem domain (0, 1] when the rotation velocity in an 
infinitesmal neighborhood around r = is modified such that V(0) becomes zero, which corresponds 
to replacing p(0) = oo with a finite (large) value of p(0). Such a slight modification of rotation curve 
results in no practical difference in the computed mass density distribution and the value of total galactic 
mass Aid, while providing great convenience for numerical computation. 

Therefore, we take the same approach here to slightly modify the rotation-curve data files of Sofue 
such that the first point at r = has V(0) = while leaving all the rest data points unchanged; the 
resulting rotation curves are shown in figure 3 with r measured in units of R g = 20.55 and 20.00 
(kpc), rotation velocity V(r) in units of Vb = 220 and 180 (km/s), respectively for the Milky Way and 
NGC 4945. Also shown in figure 3 are the computed mass density distributions for the Milky Way and 
NGC 4945. With the computed values of A = 1.6365 and 1.6873, we have M d = 1.4138 x 10 11 and 
8.9337 x 10 10 (solar-mass) respectively for the Milky Way and NGC 4965, according to (O. For the 
Milky Way, one unit of nondimensional p corresponds to as surface mass density of Md hj R 2 g — 3.35 
(solar-mass/pc 2 ). In the solar neighborhood around 8 (kpc) from the galactic center, which corresponds 
to r = 0.3893, we have p ~ 43 (from figure 3) and therefore the surface mass density around the Sun 
should be ~ 144 (solar-mass/pc 2 ). 

Due to the large central peaks in rotation curves near r = 0, the computed mass densit y profiles show 



sharp increase of p toward the galactic center, as consistent with the previous findings of iFeng & Gallo 



(201 1) based on a series of idealized rotation curves. Because the rotation curves of the Milky Way and 
NGC 4945 are generally flat, the mass density profiles show only about one order of magnitude decrease 
in the large interval (0.1, 0.9), unlike that for NGC 4736 with more than two order of magnitude decrease 
corresponding to a rotation curve of velocity generally decreasing with radial distance. 

3.3 NGC 224, NGC 5055 

The g alaxies NGC 224 and NGC 5055 were classified as those with rotation curves having "no central 

peak' ' dSofue etal.1 1999b . in contrast to that of the Milky Way. Their rotation curves (from http://www.ioa.s. u-tokyo.ac.jp/~sofue) 



and our computed mass density profiles are presented in figure 4, with r measured in units of R g = 31.25 
and 39.35 (kpc), rotation velocity V(r) in units of Vq = 250 and 190 (km/s), respectively for NGC 224 
and NGC 5055. Corresponding to the rotation curves without the central peak, the mass density profiles 
of NGC 224 and NGC 5055 vary less dramatically as r — > than those with large central peaks in 
figure 3. With the computed values of A = 1.6450 and 1.6888, we obtain M d = 2.7619 x 10 11 and 
1.9567 x 10 11 (solar-mass) respectively for NGC 224 and NGC 5055. 

3.4 NGC 2403, NGC 3198 



With their rotation curves being classified as "rigid-body type" dSofue et al.l l 1999). the NGC 2403 and 
NGC 3198 galaxies have rotation velocity increase gradually from the galactic center almost like rigid- 
body rotation for a considerable radial distance before leveling off. Figure 5 shows the rotation curves 
(from http://www.ioa.s. u-tokyo.ac.jp/~sofue) and corresponding mass density profiles of NGC 2403 
(thick line, which also has the nonzero velocity at r = replaced with V(0) = 0) and NGC 3198 
(thin line), with r measured in units of R g = 19.70 and 31.05 (kpc), rotation velocity V(r) in units of 
Vb = 130 and 160 (km/s), respectively for NGC 2403 and NGC 3198. The peak density at r = in 
figure 5 is further reduced from that in figure 4, due to less steep change in the rotation velocity around 
the galactic core. With the computed values of A — 1.4918 and 1.6022, we obtain M d = 5.1915 x 10 10 
(solar-mass) for NGC 2403 and 1.1541 x 10 11 (solar-mass) for NGC 3198. 

It is noteworthy that the NGC 3198 rotation curve has a small spike near r = 0, which results 
in a sharp turn in the mass density around the same location. Another obvious wiggling spike in the 
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r 

Figure 4: Profiles of NGC 224 and NGC 5055 rotation velocity V(r) and mass density p(r), with the thick 
line for that of NGC 224 and the thin line for NGC 5055. The computed values of the galactic rotation 
number A are 1.6450 and 1.6888 for cases of NGC 224 and NGC 5055, respectively. 
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Figure 5: Profiles of NGC 2403 and NGC 3198 rotation velocity V(r) and mass density p(r), with the thick 
line for that of NGC 2403 and the thin line for NGC 3198. The computed values of the galactic rotation 
number A are 1.4918 and 1.6022 for cases of NGC 2403 and NGC 3198, respectively. 
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rotation curve is at r ~ 0.2 causing a corresponding corner formed in the mass density profile in that 
neighborhood. Apparently the effects of some of the fine features in the rotation curve are confined 
locally in a small nearby neighborhood. 

4 Discussion 

4.1 Nonzero rotation velocity at r = 

In §3.2 we have treated the rotation curves having nonzero velocity at r = by replacing the value of 
V(0) with a zero value at this one data point in the rotation curve data file from measurements. Such a 
simplistic approach may be a little distasteful to some people of rigorous mind. So, a more elaborated 
treatment is provided here. 

With the thin-disk model, we have illustrated with galaxies of various types of realistic rotation 
curves that the mass density is always highest at the galactic center, and a nonzero rotation velocity at 
r = corresponds to an infinite mass density at the galactic center. To enable numerical treatment of 
the infinite local mass density, it may not be unreasonable to consider the galaxies with nonzero rotation 
velocity at r = to consist of a dense spherical core at the galactic center in addition to a self-gravitating 
thin disk. In that case, we should modify © to include a term due to the dense core with a spherically 
symmetric gravitational field. Among many choices, we can simply assume a spherical core confined 
within a small volume, e.g., in r < R c = 0.0001, having a mass M(r) — AV(0) 2 r where V(0) is 
nonzero according to the measured rotation curve. This corresponds to a spherically symmetric mass 
density p(r) = [dM (r) / dr}/ (2Trr 2 ) — A V (0) 2 / {2irr 2 ) in r < R c , becoming infinite as r — > 0. As 
a consequence, the second term in {§), namely V(r) 2 , can be replaced by ^A [V(r) 2 ~ V^(0) 2 ] for 
r < R c and by \A |V(r) 2 — V(0) 2 R c /r] for r > R c . Such a modification is actually equivalent to 
replacing the original rotation curve V(r) with a modified one that becomes zero at r = as 



If we apply this approach to the Milky Way, which has V(0) = 0.9282, with R c = 10~ 4 we 
obtain A = 1.6368 (instead of 1.6365 in §3.2). Thus, the total mass in the galactic disk is Md = 
1.4135 x 10 11 (solar-mass) (instead of 1.4138 x 10 11 solar-mass in §3.2). The mass in the spherical core 
is AV(0) 2 R c M d = 1.6368 x 0.9282 2 x 1.4135 x 10 7 = 1.9933 x 10 7 (solar-mass). The combined 
mass of the core and disk is then 1.4137 x 10 11 (solar-mass) that is of no practical difference from the 
value in §3.2. With such a small core of R c = 10~ 4 , the modified rotation curve (fT6] l is also not of 
practical difference from that (thick line) in figure 3, having < 2% change at r = 0.0024 (the second 
data point in measured rotation curve), < 1% change at r = 0.0049 (the third data point), ~ 0.5% 
change at r = 0.0073 (the forth data point), and so on so forth. 

However, if we take R c = 0.01 for a bigger core, the computed value of A for the Milky Way 
becomes 1.6599 and the corresponding mass in the galatic disk then is Md = 1.3939 x 10 11 (solar- 
mass). Combining with the mass of the core (A V(0) 2 R c Md = 1.9934 x 10 9 solar-mass), we have a 
total galactic mass of 1.4138 x 10 11 (solar-mass), which is basically the same as that in §3.2. Hence, 
the total galactic mass remains unchanged for a substantial range of the spherical core size R c . But 
with R c = 0.01, the modified Miky Way disk rotation curve according to (fTSl differ noticeably from 
the original one provided by Sofue (as shown in figure 6), especially around the galactic core where the 
influence of gravitational field of the spherical core is more significant. Yet the computed p(r) in the 
thin disk still appears indistinguishable from that in figure 3, except the peak value at r = is reduced 
to 3650 from 25262 (in figure 3). This is because in a small core at the galactic center the details of mass 
distribution, whether axisymmetrically or spherical symmetrically, cannot make much a difference in the 
gravitational field some distance away. 

It seems though that the effort of decomposing the galaxy into a small spherical core and a thin disk 
only helps treat rotation curves with nonzero velocity at r = in a mathematically somewhat convenient 
manner, such that the need of explicitly considering the infinite mass density is eliminated. Physically, 
a nonzero rotation velocity at r = has unclear meaning and should remain as a subject of debate; and 



{ 



y/V(r) 2 -V{0) 2 , r<R c 
y/V{r) 2 -V(0) 2 R c /r, r > R c 



(16) 
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Figure 6: Profiles of the Milky Way rotation curves when decomposed into that corresponding to a spherical 
core of R c = 0.01 and a thin disk. The original rotation curve from Sofue is also shown here as a reference. 



so should the meaning of the corresponding infinite mass density, because the common wisdom usually 
indicates that "nature abhors infinities". Thus, we prefer the straightforward treatment in §3.2 to simply 
bring the rotation velocity to zero at r = 0, especially when it does not seems to be at the expense of 
general result accuracy. The insensitivity of mass distribution in the galactic disk and total galactic mass 
to detailed descriptions of the structure in a small spherical central core illustrated here is consistent with 
the findings of some previous authors (cf . iNordsieckll 1 973l and citations therein). 



4.2 Central bulge in disk galaxy 

Yet our methodology for treating a central spherical core can be easily extended for analyzing galaxies 
with considerably larger central bulge with a priori given spherical mass distributions. As an example, 
assuming the Milky Way rotation curve in figure 3 to be a result of the combination of a central bulge 
with a spherically symmetric mass density 

Pb(r) = P me- {r/R ^ , (17) 

and an axisymmetrically distributed mass in a thin-disk, the disk portion of mass distribution then must 
satisfy (O with V(r) in figure 3 being replaced by 

\jv{rf - [1 - e-(-/*>»] . (18) 

For Rb = 0.2 and pbo/A = 7, the disk rotation curve as determined from (Tf8l is shown in figure 
7 together with the computed disk mass density distribution. In the presence of this bulge, the value 
of A becomes 2.4793 and the disk mass density shows a dip around r = 0.12. Thus, p^o = 7 x 
2.4793 = 17.3551 and the corresponding bulge density profile is also shown in figure 7. The mass 
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Figure 7: Profiles of the Milky Way rotation velocity and mass density for the disk portion and bulge portion 
(with Rh = 0..2 and pw/A = 7) as noted along with those in figure 3 (thick line) as references. 
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in the disk portion (calculated from (fT2] >) is Md = 9.3320 x 10 10 (solar-mass) and that in the bulge 
portion M b = 4n p b0 i^M d /3 = 5.4273 x 10 10 (solar-mass). The total galactic mass M g = M d + M b 
= 1.4759 x 10 11 (solar-mass) (instead of 1.4138 x 10 11 solar-mass predicted by a pure disk model in 
§3.2). Because substantial amount of the mass is concentrated in the central bulge with its portion of 
spherically symmetric mass density practically diminshes for r > 0.3, the disk surface mass density 
in the solar neighborhood around r = 0.3893 (corresponding to 8 kpc) becomes p(0.3893) Mdh/R 2 g 
= 108 (solar-mass/pc 2 ) where p(0.3893) ~ 49 from the "disk" mass density curve in figure 7. Even 
though the presence of our example bulge causes only a few percent of increase in the total galactic mass 
from that predicted by a pure disk model, the disk surface mass density in the solar neighborhood can 
decrease by 25%. 

If R b — 0.25 and pbo/A = 5, the value of A will become 3.0223 and therefore p b o = 15.1115. 
As a consequence, Md — 7.6554 x 10 10 (solar-mass), M b = 7.5716 x 10 10 (solar-mass), and M g = 
1.5227 x 10 11 (solar-mass). The value of p(0.3893) is ~ 41 corresponding to the disk surface mass 
density in the solar neighborhood of 74 (solar-mass/pc 2 ). Thus, for a given rotation curve, the actual 
value of disk surface mass density in a galaxy can vary in a wide range when considering a model 
with a combination of a spherical bulge and an axisymmetric thin disk, depending upon the bulge mass 
structure. The bulge mass structure described by (TlTt is considered here only for illustrative purpose with 
the convenience for mathematical manipulation. The fact that adding a spherical bulge offers a another 
degree of freedom to adjust mass distribution in the galactic disk should be independent of the bulge 
mass structure. This extra degree of freedom comes at the expense of uncertainty due to the difficulties 
in determining the bulge mass structure that is governed by much more complicated physical processes 
than simply balancing the gravitational force and centrifugal force. Hence we chooce to take the bulge 
mass structure as given a priori in analysing mass distribution in disk galaxies according to Newtonian 
dynamics, with our focus kept on the thin-disk portion of galaxies. 

Without considering the central bulge, the mass distribution in galactic disk can be uniquely deter- 
mined corresponding to a given rotation curve. With the central bulge, its mass structure must be known 
a priori in order to compute a unique disk mass distribution corresponding to the given rotation curve. 
But how to reliably determine the bulge mass structure besides using its luminosity information and an 
assumed mass-to-light ratio seems to still be an open question. 

Nevertheless, our illustrative analysis presented here demonstrates the general effect of a central 
bulge as to basically shift mass from periphery toward the center of a galaxy for a given rotation curve. 
The more massive a central bulge becomes, the less mass is needed in the disk periphery region according 
to Newtonian dynamics. Yet the total mass in a galaxy seems to be much less sensitive to the presence 
of a central bulge. 

4.3 Rotation velocity beyond the cut-off radius 

Also as shown in §3. 1, our finite-disk galaxy model and the associated computational method can further 
be used to determine the rotation velocity of matters outside the cut-off radius, which we assume to be 
the edge of the galaxy where the mass density diminishes. Again taking the Milky Way as an example, 
figure 8 shows the computed rotation velocity beyond the galactic edge r = 1, as a continuation from 
the measured rotation curve that ends at r = 1 and gradually approaching the Keplerian rotation curve 
for r > 2. Here the Keplerian rotation curve is generated by applying Keplerian dynamics 



with p(f) and A being obtained through computations in §3.2. Because Keplerian dynamics cannot 
correctly describe the situation of disk galaxies with nonspherically symmetic gravitational field, the ro- 
tation curve predicted by Keplerian dynamics ([T9i l from the disk mass distribution p(f) differs noticeably 
from that of Newtonian dynamics (as depicted with the thick line and its extension in figure 8). Only at 
a large distance (e.g., r > 2) from the galactic disk does the Keplerian rotation curve approaches that 
computed based on Newtonian dynamics, for the effect of disk structure diminishes at large distance 
where the gravitation field of a finite disk galaxy approaches that of a point mass. 




(19) 
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Figure 8: Profiles of the Milky Way rotation curves of the original measured one (with imposed 1^(0) = 0, 
thick line) and its extension outside the disk edge at r = 1 (thin line as labeled), as well as Vk{t) according 
to Keplerian formula ( fT9l ) based on the mass density p{r) shown in figure 3 (thin line as labeled). 
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Figure 9: Profiles of the Milky Way mass within radial distance r (thick line), M(r) = 
and that estimated with Keplerian dynamics (thin line as labeled), Mx(r) = ArV(r) 2 . 



2ir h Jq p(f) rdr, 



4.4 Applicability of Keplerian dynamics 

If Keplerian dynamics were applied to estimate the amount of mass within the solar radius (8 kpc cor- 
responding to r = 0.3893) from the measured local rotation velocity (201.0658 km/s corresponding to 
^(0.3893) = 0.9139), we would obtain M K {r) = ArV(r) 2 = 1.6365 x 0.3893 x 0.9139 2 = 0.5321 
which corresponds to 0.5321 x 1.4138 x 10 11 = 0.7523 x 10 11 (solar-mass). The actual amount of mass 
within the solar radius calculated using M (r) = 2n h J Q p(f) fdf at r = 0.3893 based on the computed 
p(r) in §3.2 is 0.5078, which corresponds to 0.7179 x 10 11 (solar-mass). Although the mass within the 
solar radius (r = 0.3893) estimated with Keplerian dynamics (Mr- (0.3893) = 0.5321) does not seem 
too far off the actual value (M (0.3893) = 0.5078), the value of Mjc(r) deviates more and more from 
M(r) with increasing r as can be seen in figure 9. The value of Mjf(r) may even decrease with r, when 
calculated according to measured rotation curve (as clearly shown in figure 8 for r > 0.85). Because 
Mk{t) is expected to monotonically increase with r for there is no physical evidence of negative mass in 
the universe, a negative slope of Mr- if) versus r indicates a failure of Keplerian dynamics for correctly 
deriving the mass distribution corresponding to measured rotation curve for the Milky Way. Or, in other 
words, a rotation curve that does not satisfy d[r V(r) 2 ]/dr > 0, namely 



dV(r) 
dr 



> 



V(r) 
2r 



(20) 



is inconsistent with spherically symmetric gravitational potential and thus Keplerian dynamics becomes 
inapplicable in a strict sense. 



The same condition of d20Ti was referred t o as the spheric i ty con dition bv lJalocha. Brateck & Kutschera 



( 2008 ): Bratek. Jalocha & Kutschera ( 2008 ): Jalocha. et al. ( 2010l) . and the violation of which was used 
as an indication of disk model being more appropriate for determining the mass distribution and the 
presence of a massive spherical halo of non-baryonic dark matter being unlikely. Actually ( |20i > is only 
a necessary condition for the sphericity of gravitational potential to exist, but not sufficient. A rotation 
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curve satisfying (|20T> does n ot guarantee that it must correspond to a spherically symmetric gravitational 



field. Feng & Galloj d201 II) showed that a flat rotation curve can be described by both a spherically 



symmetric and an axisymmetric disk mass distribution. But using a spherically symmetric mass model, 
namely Keplerian dynamics, to describe a rotating disk galaxy can lead to erroneous results and conclu- 



4.5 Circular orbit stability 

It is interesting to note that the mathematical form of the sphericity condition (|20| > appears quite similar 
to the necessary condition for circular orbit stability 

«m > - m , (2D 

dr r 

as can be derived from the consideration of angular momentum conservation for a rotating object slightly 
deviating from its original (circular) orbit as follows. An object that is rotating with a velocity V(r) at 
radial coordinate r possesses an angular momentum r V(r). If it deviates its original orbit at r to r + Sr 
(due to some sort of perturbations), its rotation velocity should change from V to V + SV such that 
(r + Sr) {V + SV) = r V or SV/Sr = -V/r (< 0, i.e., 5V < when Sr > and SV > when 
Sr < 0) according to the conservation of angular momentum. On the other hand, this object is subjected 
to a gravitational force at r + Sr equal to mV(r + Sr) 2 / (r + Sr), where m is its mass and V(r + Sr) 
is the rotating velocity of objects at r + Sr according to the rotation curve. Thus, for this object to be 
pulled back by gravitational force to its original orbit, namely for its orbit to be centrifugally stable, we 
must have V + SV < V(r + Sr) for Sr > and V + SV > V(r + Sr) for Sr < 0. This leads to SV 
< Sr dV(r)/dr + o(Sr 2 ) for Sr > and SV > Sr dV(r)/dr + o{Sr 2 ) for Sr < as a result of Taylor 
expansion around Sr = 0, namely SV/Sr = —V/r < dV/dr as Sr — > 0. 

In the case of the solar system with a point mass at r = 0, the planets rotate following the Keplerian 
rotation curve dV / dr = —V/ (2r) (taking the equal sign in d20l > for the mass Mjc(r) does not change for 
r > 0). Because —1/2 > —1, the Keplerian rotation curve satisfies the circular orbit stability condition 
(|2TT > as evidenced by the existence of the solar system with many planets circling around the Sun year 

after year. 

Many spiral galaxies exhibit nearly flat rotation curves (cf . ISofue & Rubnl boO 1 ) corresponding to 



dV/dr ~ which can easily satisfy (f2TT> . Thus, the rotating matters distributed in circular orbits of 
galactic disk, as can be computed with the method illustrated in the present work, for flat rotation curves 
are stable in the sense of that similar to the planets circling around the Sun. But d2"TT i is only a necessary 
condition for stability. There have been many other (necessary) conditions proposed i n the literature 



for rot ating disk galaxy stability, which often seem controversial as critically discussed bv ljalocha. et al 



( 20ld) . The circular orbit stability condition (OTt derived in the present work is established from concrete 



physical principle and can be used to examine the validity of measured rotation curves. Especially for 
those rotation curves containing decreasing velocity at large radius, the stability condition (I2H with 
its right side oc — 1/r is likely violated. The portion of a rotation curve not satisfying d2"TT i may point 
to the issues with too serious deviations from circular orbits and axisymmetry due to the spiral arms. 
After all, the axisymmetric disk model with rotation velocity depending only on radius is a tremendous 
simplification of a realistic rotating galaxy; such simplified description of the reality should be constantly 
checked for consistency. 

For the Milky Way rotation curve (cf. figure 3), the portion of r > 0.8 has dV/dr ~ —0.75 while 
0.9 < V/r < 1.32. Thus, we have —0.75 > —0.9 that satisfies the circular orbit stability condition (|2H 
but not the sphericity condition (l20b which consistent with that shown in figure 9. Thus, the Milky Way 
appears to be appropriately described with the thin-disk model or with a combination of a central bulge 
and a thin disk. 

For the NGC 4736 rotation curve of Sofue (cf. figure 1), the negative slope in the outer region 
r > 0.6 is dV/dr ~ —0.875 while 0.85 < V/r < 2. Thus, circular orbit instability is likely to occur 
in the outer region of NGC 4736 if the rotating matters indeed foll ow the rotation curv e of Sofue. In 
a relative sense, the THINGS version of NGC 4736 rotation curve Jde Blok et al.ll2008h has less steep 



negative slope than that of Sofue, indicating that the THINGS version describes more stable circular 
motion of rotating matters. 
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5 Conclusions 



With the computational method presented here, the mass distributions in mature spiral galaxies corre- 
sponding to various types of measured rotation curves can be accurately determined, based on Newtonian 
dynamics without the need of fictitious rotation velocity outside the cut-off radius. In our finite thin-disk 
model, the galactic disk edge where the mass density precipitously diminishing is assumed to coincide 
with the cut-off radius in rotation curve measurement, based on the physical intuition that inability of 
measurement is a co nsequence of the absence of matter outside the cut-off radius. In view of some 
authors' desires (e.g.. lJalocha. Brateck & Kutscherall2008i in the case of NGC 4736) in considering the 
small amount of hydrogen HI detected outside the cut-off radius, we also examine the effect of such 
HI with an a posteriori computation by extending our computational domain beyond the disk edge. Not 
surprisingly, our result shows that the HI outside the cut-off radius only has inconsequential effect on the 
rotation velocity because it accounts for only a fraction of a percent of the mass in the galaxy. Therefore, 
our intuitive assumption of galactic disk edge at the cut-off radius is not expected to result in numerical 
errors of any practical significance, unless the measured cut-off radius is too far off the true edge of the 
galaxy. 

Despite the difficulities in clarifying the physical meaning, nonzero rotation velocities at the g alactic 
center r = were reported in rotation curve measurements for several galaxies dSofue et al . 1999). The 
nonzero value of rotation velocity at r = mathematically corresponds to unbounded local mass density, 
namely p(0) oo which is intractable in numerical computations. Our first choice of the possible 
approaches is to simply replace the nonzero rotation velocity with 1^(0) for the measurement data point 
at r = in pure thin-disk model and compute the corresponding mass distribution in the thin-disk 
galaxy. Our second choice is to place a small spherical core at r = that avoids explicit mathematical 
consideration of the infinite local mass density. Because the mass in spherical core contributes to the 
gravitation force in the galactic disk, the rotation velocity in galactic disk is modified accordingly by 
subtracting out the spherical core effect and the disk mass distribution can be computed based on the 
modified rotation curve. Our results show that as long as the spherical core is small (e.g., being confined 
within r = 0.01), no noticeable change can be observed in the general disk mass distribution except the 
nonzero rotation velocity at r = can be mathematically dealt with. 

Observations has shown that many mature spiral galaxies have bright central bulges of various sizes 
in addition to their disk-like mass distribution. To examine the basic effects of a central bulge, we assume 
the bulge has a spherically symmetric mass structure such that its gravitational effect can be conveniently 
incorporated in our thin-disk model formulation with a modified rotation curve for the disk portion. The 
computational results suggest that the presence of a central bulge tends to effectively shift the mass from 
periphery toward the galactic center with little change in the total galactic mass. 

Extending the computational domain beyond the galactic edge enables us to also compute rotation 
velocity outside the cut-off radius. Outside the galactic edge where we assume the amount of mass is 
negligible, the computed rotation velocity does not exactly follow the Keplerian rotation velocity until 
out over r > 2. With the mass distributed in a thin disk as determined according to Newtonina dynamics, 
the computed Keplerian rotation velocity within the galactic disk differs noticeably from that of measured 
rotation curve. 

By applying the principle of angular momentum conservation, we can derive a necessary condition 
for circular orbit stability. It appears that the galaxies with flat or increasing rotation velocities with radius 
are more stable due to angular momentum conservation than those with decreasing rotation velocities. 
Especially in the region near the galactic edge, the rotation curves having too steep of negative slope may 
violate the condition for circular orbit stability and therefore their validity for realistically describing the 
galactic rotational characteristics may become questionable. 
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A Treatments of Singular Elements 

As in lFeng&Gallol(l20l 3), the complete elliptic integrals of the first k ind and second kind in (0 can be 
numerically computed with the formulas (lAbramowitz &Stegunll 19721) 



K(m) — aim\ — log(rai) bim\ 



1=0 



and 



where 



E(m) = 1 + ^2 Cim i ~ lo s( TO i) X! dim i ' 



(22) 
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(24) 
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Clearly, the terms associated with K(m,i) and E(rrii) in © become singular when f 
elements with 7"j as one of their end points. 

The logarithmic singularity can be treated by converting the singular one-dimensional integrals into 
non-singular two-dimensional integrals by virtue of the identities: 



Jo /(£) log^£ = - Jo Jo m)dr)<% 
So f(0 log(l - 0« = - Jo Jo /(! - Wdvdt 



(25) 



where /(£) denotes a well-behaving (non-singular) function of £ on < £ < 1. 

However, a more serious non-integrable singularity l/(f — r,-) exists due to the term E(rrii)/ (f — r,) 
in (O as f — » r^. The l/(f — r, ) type of singularity is treated by taking the Cauchy principle value to 
obtain meaningful evaluation fcf.|Kanwa]|ll99a). as commonly done with the boundary element method 
dSladek & Sladek 1998; Sutradhar, Paulino & Grav 2008). In view of the fact that each r; is considered 
to be shared by two adjacent elements covering the intervals [rj_i, ri\ and [rj, the Cauchy principle 
value of the integral over these two elements is given by 



lim 

e-s-0 



p(f)fdf 



(26) 



In terms of elemental £, d26] l is equivalent to 



lim 

£^0 



1-C 

[Ml - £) + - + r- i+ i£]d£ ' 



e/(r <+I -r<) 

Performing integration by parts on (l27t yields 



£ 



(27) 
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where the two terms associated with log e cancel out each other, the terms with e log e become zero at the 
limit of e — 5- 0, and the first term becomes nonzero when the mesh nodes are not uniformly distributed 
(namely, the adjacent elements are not of the same segment size). 
At the galaxy center r, = 0, 



Thus, the — r<) type of singularity disappears naturally. However, num erical difficul ty can still arise 
if p itself becomes singular as r — > 0, e.g., p oc 1/r as for the Mestel disk dMestellll963b . The singular 
mass density at r = corresponds to a mathematical cusp, which usually indicates the need of finer 
resolution in the physical space. To avoid the cusp in mass density at the galactic center, we can impose 
a requirement of continuity of the derivative of p at the galaxy center r = 0. This be easily implemented 
at the first node i = 1 to demand dp/dr = at r = 0. In discretized form for r% = we simply have 



When n = 1 (i.e., i = N), we are at the end node of the problem domain. Here we use a numerically 
relaxing boundary condition by considering an additional element beyond the domain boundary covering 
the interval [n, r<+i], because it is needed to obtain a meaningful Cauchy principle value. In doing so we 
can also assume r^+i— n = r,— Ti-% such that log[(rj+i — r i)/( r i — r i-i)\ becomes zero, to simplify the 
numerical implementation. Moreover, it is reasonable to assume pi + i = because it is located outside 
the disk edge. With sufficiently fine local discretization, this extra element can be considered to cover a 
diminishing physical space such that its existence becomes numerically inconsequential. Thus, at = 1 
we have 



Now that only logarithmic singularities are left, (l25t can be used to eliminate all singularities in comput- 
ing the integrals in (|5). 
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